convert p values to qvalues

#read in all phase1 significance values:
all.sigs<-read.csv(file="~/Desktop/anoph.3.march.2020/all.sigs.csv")
table(all.sigs$chrom)
## 
##     2L     2R     3L     3R      X 
## 638059 946059 524268 737401 198852
#convert 0 bayescenv q-vals to 1e-4, to avoid infinite values when plotting with -log10()
all.sigs$prec.q[all.sigs$prec.q == 0]<-.0001
all.sigs$temp.q[all.sigs$temp.q == 0]<-.0001
all.sigs$hum.q[all.sigs$hum.q == 0]<-.0001

#convert lfmm prec p values to genome-wide FDR adjusted q values
for (i in levels(all.sigs$chrom)){
qvals<-qvalue(all.sigs$prec.p[all.sigs$chrom== i])
all.sigs$prec.p[all.sigs$chrom== i]<-qvals[["qvalues"]]
}
#convert lfmm hum p values to genome-wide FDR adjusted q values
for (i in levels(all.sigs$chrom)){
  qvals<-qvalue(all.sigs$hum.p[all.sigs$chrom== i])
  all.sigs$hum.p[all.sigs$chrom== i]<-qvals[["qvalues"]]
}
#convert lfmm temp p values to genome-wide FDR adjusted q values
for (i in levels(all.sigs$chrom)){
  qvals<-qvalue(all.sigs$temp.p[all.sigs$chrom== i])
  all.sigs$temp.p[all.sigs$chrom== i]<-qvals[["qvalues"]]
}
rm(qvals) #space saver

#give order for plotting manhattan plots
all.sigs$order<-c(all.sigs$pos[all.sigs$chrom== "2R"],
                  all.sigs$pos[all.sigs$chrom== "2L"]+61545105,
                  all.sigs$pos[all.sigs$chrom== "3R"]+61545105+49364325,
                  all.sigs$pos[all.sigs$chrom== "3L"]+61545105+49364325+53200684,
                  all.sigs$pos[all.sigs$chrom== "X"]+61545105+49364325+53200684+41963435)

lfmm manhattans

#LFMM precip plot -log10(signficance vals) as manhattan plots
plot1<-ggplot(data=all.sigs, aes(x=order, y=-log10(prec.p), col = chrom))+
  geom_scattermore(pointsize = 1.2)+
  theme(panel.grid=element_blank(), panel.background=element_blank())+
  geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant precip outliers
lfmm.prec.outliers<- all.sigs[all.sigs$prec.p < .01,]
nrow(lfmm.prec.outliers) #check number of outliers
## [1] 17113
#LFMM hum plot -log10(signficance vals) as manhattan plots
plot2<-ggplot(data=all.sigs, aes(x=order, y=-log10(hum.p), col = chrom))+
  geom_scattermore(pointsize = 1.2)+
  theme(panel.grid=element_blank(), panel.background=element_blank())+
  geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant precip outliers
lfmm.hum.outliers<- all.sigs[all.sigs$hum.p < .01,]
nrow(lfmm.hum.outliers) #check number of outliers
## [1] 10981
#LFMM temp plot -log10(signficance vals) as manhattan plots
plot3<-ggplot(data=all.sigs, aes(x=order, y=-log10(temp.p), col = chrom))+
  geom_scattermore(pointsize = 1.2)+
  theme(panel.grid=element_blank(), panel.background=element_blank())+
  geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant precip outliers
lfmm.temp.outliers<- all.sigs[all.sigs$temp.p < .01,]
nrow(lfmm.temp.outliers) #check number of outliers
## [1] 12396
#
grid.arrange(plot1,plot2,plot3, nrow=3)

bayescenv manhattans

#bayescenv precip plot -log10(signficance vals) as manhattan plots
plot1<-ggplot(data=all.sigs, aes(x=order, y=-log10(prec.q), col = chrom))+
  geom_scattermore(pointsize = 1.2)+
  theme(panel.grid=element_blank(), panel.background=element_blank())+
  geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant precip outliers
baye.prec.outliers<- all.sigs[all.sigs$prec.q < .01,]
nrow(baye.prec.outliers) #check number of outliers
## [1] 3857
#bayescenv hum plot -log10(signficance vals) as manhattan plots
plot2<-ggplot(data=all.sigs, aes(x=order, y=-log10(hum.q), col = chrom))+
  geom_scattermore(pointsize = 1.2)+
  theme(panel.grid=element_blank(), panel.background=element_blank())+
  geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant hum outliers
baye.hum.outliers<- all.sigs[all.sigs$hum.q < .01,]
nrow(baye.hum.outliers) #check number of outliers
## [1] 3688
#bayescenv temp plot -log10(signficance vals) as manhattan plots
plot3<-ggplot(data=all.sigs, aes(x=order, y=-log10(temp.q), col = chrom))+
  geom_scattermore(pointsize = 1.2)+
  theme(panel.grid=element_blank(), panel.background=element_blank())+
  geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant temp outliers
baye.temp.outliers<- all.sigs[all.sigs$temp.q < .01,]
nrow(baye.temp.outliers) #check number of outliers
## [1] 10195
grid.arrange(plot1,plot2,plot3, nrow=3)

par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(na.omit(lfmm.hum.corrs$id[abs(lfmm.hum.corrs$diff) < .05]))
sig.pos<-sig.pos[1:27]
for (i in sig.pos){
  frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
  frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
  plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
  abline(lm(frq~sort.pops.phase1$hannual))
  points(sort.pops.phase2$hannual,frq2, col="red")
  abline(lm(frq2~sort.pops.phase2$hannual), col="red")
}

par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(na.omit(lfmm.hum.corrs$id[abs(lfmm.hum.corrs$diff) > .8]))
sig.pos<-sig.pos[1:27]
for (i in sig.pos){
  frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
  frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
  plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
  abline(lm(frq~sort.pops.phase1$hannual))
  points(sort.pops.phase2$hannual,frq2, col="red")
  abline(lm(frq2~sort.pops.phase2$hannual), col="red")
}

check out the correlations of allele frequency patterns with humidity

#visualize
plot(lfmm.hum.corrs$diff)
abline(h=0, col="red")

plot(lfmm.hum.corrs$phase1, lfmm.hum.corrs$phase2)

hist(lfmm.hum.corrs$diff)

plot(-log10(lfmm.hum.outliers$hum.p), all.cases$diff)

plot(-log10(lfmm.hum.outliers$hum.q), all.cases$diff)

plot(lfmm.hum.outliers$hum.fst, all.cases$diff)

#different classes of relationships?

layout(matrix(c(1,1,1,1,2,3,4,5),2, 4, byrow=TRUE))
plot(lfmm.hum.corrs$phase1)
abline(h=c(.4,.1,-.25), col="red")
hist(lfmm.hum.corrs$diff[lfmm.hum.corrs$phase1 > .4])
abline(v=c(-.2,.2), col="red")
hist(lfmm.hum.corrs$diff[lfmm.hum.corrs$phase1 < .4 & lfmm.hum.corrs$phase1 > .1])
abline(v=c(-.2,.2), col="red")
hist(lfmm.hum.corrs$diff[lfmm.hum.corrs$phase1 < .1 & lfmm.hum.corrs$phase1 > -.25])
abline(v=c(-.2,.2), col="red")
hist(lfmm.hum.corrs$diff[lfmm.hum.corrs$phase1 < -.25])
abline(v=c(-.2,.2), col="red")

#strongly positive correlation

par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(lfmm.hum.corrs$id[lfmm.hum.corrs$phase1 > .4])
sig.pos<-sample(sig.pos, size=27, replace = FALSE)
for (i in sig.pos){
  frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
  frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
  plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
  abline(lm(frq~sort.pops.phase1$hannual))
  points(sort.pops.phase2$hannual,frq2, col="red")
  abline(lm(frq2~sort.pops.phase2$hannual), col="red")
}

#weakly positive correlation

par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(lfmm.hum.corrs$id[lfmm.hum.corrs$phase1 < .4 & lfmm.hum.corrs$phase1 > .1])
sig.pos<-sample(sig.pos, size=27, replace = FALSE)
for (i in sig.pos){
  frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
  frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
  plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
  abline(lm(frq~sort.pops.phase1$hannual))
  points(sort.pops.phase2$hannual,frq2, col="red")
  abline(lm(frq2~sort.pops.phase2$hannual), col="red")
}

#weakly negative correlation

par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(lfmm.hum.corrs$id[lfmm.hum.corrs$phase1 < .1 & lfmm.hum.corrs$phase1 > -.25])
sig.pos<-sample(sig.pos, size=27, replace = FALSE)
for (i in sig.pos){
  frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
  frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
  plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
  abline(lm(frq~sort.pops.phase1$hannual))
  points(sort.pops.phase2$hannual,frq2, col="red")
  abline(lm(frq2~sort.pops.phase2$hannual), col="red")
}

#strongly negative correlation

par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(lfmm.hum.corrs$id[lfmm.hum.corrs$phase1 < -.25])
sig.pos<-sample(sig.pos, size=27, replace = FALSE)
for (i in sig.pos){
  frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
  frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
  plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
  abline(lm(frq~sort.pops.phase1$hannual))
  points(sort.pops.phase2$hannual,frq2, col="red")
  abline(lm(frq2~sort.pops.phase2$hannual), col="red")
}